source('../env.R')
It seems reasonable to expect that cities with simialr regional pools will have similar species entering the city, and thus a similar response to urbanisation.
city_effort = read_csv(filename(CITY_DATA_OUTPUT_DIR, 'city_effort.csv'))
Rows: 342 Columns: 7── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): city_name
dbl (6): city_id, total_city_checklists, total_city_locations, total_city_effort, total_city_area_m2, percentage_total_city_area_surveyed
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
city_effort
communities = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'communities_for_analysis.csv'))
Rows: 2462 Columns: 12── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (5): city_name, jetz_species_name, seasonal, presence, origin
dbl (4): city_id, distance_to_northern_edge_km, distance_to_southern_edge_km, relative_abundance_proxy
lgl (3): present_urban_high, present_urban_med, present_urban_low
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
communities_summary = communities %>% group_by(city_id) %>% summarise(
regional_pool_size = n(),
urban_pool_size = sum(relative_abundance_proxy > 0)
) %>% left_join(city_effort %>% dplyr::select(city_id, percentage_total_city_area_surveyed))
Joining with `by = join_by(city_id)`
ggplot(communities %>% filter(relative_abundance_proxy > 0), aes(x = relative_abundance_proxy)) + geom_bar(stat = "bin")
city_points = st_centroid(read_sf(filename(CITY_DATA_OUTPUT_DIR, 'city_selection.shp')))
Warning: st_centroid assumes attributes are constant over geometriesWarning: st_centroid does not give correct centroids for longitude/latitude data
community_data_metrics = read_csv(filename(COMMUNITY_OUTPUT_DIR, 'community_assembly_metrics_using_relative_abundance.csv')) %>%
dplyr::select(city_id, mntd_normalised, fdiv_normalised, mass_fdiv_normalised, handwing_index_fdiv_normalised, trophic_trait_fdiv_normalised, gape_width_fdiv_normalised) %>%
left_join(read_csv(filename(CITY_DATA_OUTPUT_DIR, 'realms.csv'))) %>%
left_join(communities_summary) %>%
left_join(city_points[,c('city_id', 'city_nm')]) %>%
rename(city_name='city_nm') %>%
na.omit() %>%
arrange(city_id)
Rows: 341 Columns: 43── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
dbl (43): mntd_normalised, mntd_actual, mntd_min, mntd_max, mntd_mean, mntd_sd, fdiv_normalised, fdiv_actual, fdiv_min, fdiv_max, fdiv_mean, fdiv_sd, mass_fdi...
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Rows: 342 Columns: 2── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): core_realm
dbl (1): city_id
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`Joining with `by = join_by(city_id)`
community_data_metrics
community_data_metrics %>% group_by(core_realm) %>% summarise(total_cities = n())
test_value = function(name, normalised_list) {
wilcox_test_result = wilcox.test(normalised_list, mu = 0.5)
significance = ifelse(wilcox_test_result$p.value < 0.0001, '***',
ifelse(wilcox_test_result$p.value < 0.001, '**',
ifelse(wilcox_test_result$p.value < 0.01, '*', '')))
m = mean(normalised_list)
paste(name, 'mean', round(m, 2), significance)
}
test_value('Global MNTD', community_data_metrics$mntd_normalised)
[1] "Global MNTD mean 0.5 "
test_value('Global beak gape FDiv', community_data_metrics$gape_width_fdiv_normalised)
[1] "Global beak gape FDiv mean 0.59 ***"
test_value('Global HWI FDiv', community_data_metrics$handwing_index_fdiv_normalised)
[1] "Global HWI FDiv mean 0.61 ***"
nrow(community_data_metrics)
[1] 319
Load trait data
traits = read_csv(filename(TAXONOMY_OUTPUT_DIR, 'traits_jetz.csv'))
Rows: 304 Columns: 6── Column specification ──────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────────
Delimiter: ","
chr (1): jetz_species_name
dbl (5): gape_width, trophic_trait, locomotory_trait, mass, handwing_index
ℹ Use `spec()` to retrieve the full column specification for this data.
ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(traits)
fetch_normalised_traits = function(required_species_list) {
required_traits = traits %>% filter(jetz_species_name %in% required_species_list)
required_traits$gape_width_normalised = normalise(required_traits$gape_width, min(required_traits$gape_width), max(required_traits$gape_width))
required_traits$handwing_index_fdiv_normalised = normalise(required_traits$handwing_index, min(required_traits$handwing_index), max(required_traits$handwing_index))
required_traits$mass_normalised = normalise(required_traits$mass, min(required_traits$mass), max(required_traits$mass))
traits_normalised_long = required_traits %>% pivot_longer(cols = c('gape_width_normalised', 'handwing_index_fdiv_normalised', 'mass_normalised'), names_to = 'trait', values_to = 'normalised_value') %>% dplyr::select(jetz_species_name, trait, normalised_value)
traits_normalised_long$trait = factor(traits_normalised_long$trait, levels = c('gape_width_normalised', 'handwing_index_fdiv_normalised', 'mass_normalised'), labels = c('Gape Width', 'HWI', 'Mass'))
traits_normalised_long
}
fetch_normalised_traits(c('Aplopelia_larvata', 'Chalcophaps_indica', 'Caloenas_nicobarica'))
Read in our phylogeny
phylo_tree = read.tree(filename(TAXONOMY_OUTPUT_DIR, 'phylogeny.tre'))
ggtree(phylo_tree, layout='circular')
Load resolve ecoregions
resolve = read_resolve()
Warning: st_buffer does not correctly buffer longitude/latitude datadist is assumed to be in decimal degrees (arc_degrees).
Warning: st_simplify does not correctly simplify longitude/latitude data, dTolerance needs to be in decimal degrees
to_species_matrix = function(filtered_communities) {
filtered_communities %>%
dplyr::select(city_id, jetz_species_name) %>%
distinct() %>%
mutate(present = TRUE) %>%
pivot_wider(
names_from = jetz_species_name,
values_from = "present",
values_fill = list(present = F)
) %>%
tibble::column_to_rownames(var='city_id')
}
community_nmds = function(filtered_communities) {
species_matrix = to_species_matrix(filtered_communities)
nmds <- metaMDS(species_matrix, k=2, trymax = 30)
nmds_result = data.frame(vegan::scores(nmds)$sites)
nmds_result$city_id = as.double(rownames(nmds_result))
rownames(nmds_result) = NULL
nmds_result
}
https://www.datacamp.com/tutorial/k-means-clustering-r
scree_plot = function(community_nmds_data) {
# Decide how many clusters to look at
n_clusters <- min(10, nrow(community_nmds_data) - 1)
# Initialize total within sum of squares error: wss
wss <- numeric(n_clusters)
set.seed(123)
# Look over 1 to n possible clusters
for (i in 1:n_clusters) {
# Fit the model: km.out
km.out <- kmeans(community_nmds_data[,c('NMDS1','NMDS2')], centers = i, nstart = 20)
# Save the within cluster sum of squares
wss[i] <- km.out$tot.withinss
}
# Produce a scree plot
wss_df <- tibble(clusters = 1:n_clusters, wss = wss)
scree_plot <- ggplot(wss_df, aes(x = clusters, y = wss, group = 1)) +
geom_point(size = 4) +
geom_line() +
geom_hline(linetype="dashed", color = "orange", yintercept = wss) +
scale_x_continuous(breaks = c(2, 4, 6, 8, 10)) +
xlab('Number of clusters')
scree_plot
}
cluster_cities = function(city_nmds, cities_community_data, centers) {
set.seed(123)
kmeans_clusters <- kmeans(city_nmds[,c('NMDS1', 'NMDS2')], centers = centers, nstart = 20)
city_nmds$cluster = kmeans_clusters$cluster
cities_community_data %>% left_join(city_nmds) %>% mutate(cluster = as.factor(cluster))
}
plot_nmds_clusters = function(cluster_cities) {
cluster_cities %>% dplyr::select(city_id, city_name, NMDS1, NMDS2, cluster) %>% distinct() %>%
ggplot(aes(x = NMDS1, y = NMDS2, colour = cluster)) + geom_point() + geom_label_repel(aes(label = city_name))
}
get_presence_cell_width = function(city_cluster_data_metrics) {
10 * length(unique(city_cluster_data_metrics$city_id))
}
get_presence_cell_height = function(city_cluster_data_metrics) {
species = species_in_cluster = communities %>%
filter(city_id %in% city_cluster_data_metrics$city_id) %>%
dplyr::select(jetz_species_name) %>%
distinct()
10 * nrow(species)
}
city_metric_height = 30
traits_width = 50
phylo_tree_width = 125
title_height = 8
get_image_height = function(city_cluster_data_metrics) {
get_presence_cell_height(city_cluster_data_metrics) + (3 * city_metric_height) + title_height
}
get_image_width = function(city_cluster_data_metrics) {
get_presence_cell_width(city_cluster_data_metrics) + traits_width + phylo_tree_width
}
species_in_city_label = function(species) {
paste(
ifelse(species$seasonal == 'Resident', '', substring(species$seasonal, 0, 1)),
ifelse(species$origin == 'Native', '', substring(species$origin, 0, 1)),
ifelse(species$distance_to_northern_edge_km > 200, '', paste('NRL', round(species$distance_to_northern_edge_km), sep = ' ')),
ifelse(species$distance_to_southern_edge_km > 200, '', paste('SRL', round(species$distance_to_northern_edge_km), sep = ' ')),
sep = '\n'
)
}
species_in_city_label(head(communities))
[1] "\nI\n\n" "\nI\n\n" "\nI\n\n" "P\n\n\n" "\n\n\n" "\nI\n\n"
plot_city_cluster = function(city_cluster_data_metrics, title) {
species_in_cluster = communities %>%
filter(city_id %in% city_cluster_data_metrics$city_id) %>%
dplyr::select(jetz_species_name, city_name, relative_abundance_proxy, seasonal, origin, distance_to_northern_edge_km, distance_to_southern_edge_km)
species_in_cluster$label = species_in_city_label(species_in_cluster)
tree_cropped <- ladderize(drop.tip(phylo_tree, setdiff(phylo_tree$tip.label, species_in_cluster$jetz_species_name)))
gg_tree = ggtree(tree_cropped)
gg_presence = ggplot(species_in_cluster, aes(x=city_name, y=jetz_species_name)) +
geom_tile(aes(fill=relative_abundance_proxy)) +
geom_text(aes(label=label), size=0.75) +
scale_fill_gradientn(colours=c("#98FB98", "#FFFFE0", "yellow", "orange", "#FF4500", "red", "red"), values=c(0, 0.00000000001, 0.1, 0.25, 0.5, 0.75, 1), na.value = "transparent") +
theme_minimal() + xlab(NULL) + ylab(NULL) +
theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1)) +
labs(fill='Urban Proxy Abundance')
species_in_cluster_traits = fetch_normalised_traits(species_in_cluster$jetz_species_name)
gg_traits = ggplot(species_in_cluster_traits, aes(x = trait, y = jetz_species_name, size = normalised_value)) + geom_point() + theme_minimal() + theme(axis.text.x = element_text(angle = 90, vjust = 0.5, hjust=1), axis.text.y=element_blank()) + xlab(NULL) + ylab(NULL) + labs(size = "Normalised Value")
gg_cities_mntd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = mntd_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("MNTD") + ylim(0, 1)
gg_cities_gape_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = gape_width_fdiv_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("Gape") + ylim(0, 1)
gg_cities_loco_fd = ggplot(city_cluster_data_metrics, aes(x = city_name, y = handwing_index_fdiv_normalised)) + geom_bar(stat = "identity") + theme_minimal() + theme(legend.position = "none", axis.text.x=element_blank()) + xlab(NULL) + ylab("HWI") + ylim(0, 1)
gg_title = ggplot() + labs(title = title, subtitle = paste(
test_value('MNTD', city_cluster_data_metrics$mntd_normalised),
test_value('Compound FDiv', city_cluster_data_metrics$fdiv_normalised),
test_value('HWI FDiv', city_cluster_data_metrics$handwing_index_fdiv_normalised),
test_value('Gape width FDiv', city_cluster_data_metrics$gape_width_fdiv_normalised),
sep = '\n'
)) + theme_minimal() + theme(plot.subtitle=element_text(size=8, hjust=0, color="#444444"))
gg_presence_height = get_presence_cell_height(city_cluster_data_metrics)
gg_presence_width = get_presence_cell_width(city_cluster_data_metrics)
gg_presence %>% insert_top(gg_cities_loco_fd, height = (city_metric_height / gg_presence_height)) %>% insert_top(gg_cities_gape_fd, height = (city_metric_height / gg_presence_height)) %>% insert_top(gg_cities_mntd, height = (city_metric_height / gg_presence_height)) %>% insert_left(gg_tree, width = (phylo_tree_width / gg_presence_width)) %>% insert_right(gg_traits, width = (traits_width / gg_presence_width)) %>% insert_top(gg_title, height = (title_height / gg_presence_height))
}
REGION_DEEP_DIVE_FIGURES_OUTPUT = mkdir(FIGURES_OUTPUT_DIR, 'appendix_regional_deep_dive_using_abundance')
nearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Nearctic')
nearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "San Jose" "Los Angeles" "Concord" "Tijuana" "Bakersfield"
[6] "Fresno" "Sacramento" "Mexicali" "Hermosillo" "Las Vegas"
[11] "Phoenix" "Tucson" "Durango" "Portland" "Chihuahua"
[16] "Aguascalientes" "Seattle" "Ciudad Juárez" "San Luis PotosÃ" "Mexico City"
[21] "Saltillo" "Vancouver" "Salt Lake City" "Albuquerque" "Monterrey"
[26] "Nuevo Laredo" "San Antonio" "Denver" "Austin" "Houston"
[31] "Dallas" "Oklahoma City" "Calgary" "New Orleans" "Kansas City"
[36] "Omaha" "St. Louis" "Bradenton" "Tampa" "Minneapolis [Saint Paul]"
[41] "Atlanta" "Orlando" "Louisville" "Chicago" "Indianapolis"
[46] "Milwaukee"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
nearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% nearctic_cities_community_data$city_id))
Run 0 stress 0.1005678
Run 1 stress 0.1000046
... New best solution
... Procrustes: rmse 0.007231129 max resid 0.03470438
Run 2 stress 0.1007176
Run 3 stress 0.1017503
Run 4 stress 0.1017503
Run 5 stress 0.1217441
Run 6 stress 0.1017503
Run 7 stress 0.121519
Run 8 stress 0.1210168
Run 9 stress 0.1017503
Run 10 stress 0.1000046
... New best solution
... Procrustes: rmse 0.000008641127 max resid 0.00003146381
... Similar to previous best
Run 11 stress 0.1217145
Run 12 stress 0.1212382
Run 13 stress 0.1323942
Run 14 stress 0.102946
Run 15 stress 0.1017503
Run 16 stress 0.1017503
Run 17 stress 0.1000046
... Procrustes: rmse 0.000005711568 max resid 0.00001968977
... Similar to previous best
Run 18 stress 0.1217145
Run 19 stress 0.1005678
Run 20 stress 0.1007176
*** Best solution repeated 2 times
nearctic_cities_nmds
scree_plot(nearctic_cities_nmds)
nearctic_cities = cluster_cities(city_nmds = nearctic_cities_nmds, cities_community_data = nearctic_cities_community_data, centers = 4)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(nearctic_cities)
nearctic_biomes = st_crop(resolve[resolve$REALM == 'Nearctic',c('REALM')], xmin = -220, ymin = 0, xmax = 0, ymax = 70)
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries
ggplot() +
geom_sf(data = nearctic_biomes, aes(geometry = geometry)) +
geom_sf(data = nearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_clusters.jpg'))
Saving 7.29 x 4.51 in image
test_value('Nearctic MNTD', nearctic_cities$mntd_normalised)
Warning: cannot compute exact p-value with ties
[1] "Nearctic MNTD mean 0.6 "
test_value('Nearctic beak gape FDiv', nearctic_cities$gape_width_fdiv_normalised)
Warning: cannot compute exact p-value with ties
[1] "Nearctic beak gape FDiv mean 0.59 "
test_value('Nearctic HWI FDiv', nearctic_cities$handwing_index_fdiv_normalised)
Warning: cannot compute exact p-value with ties
[1] "Nearctic HWI FDiv mean 0.31 **"
nrow(nearctic_cities)
[1] 46
nearactic_cluster1 = nearctic_cities %>% filter(cluster == 1)
plot_city_cluster(nearactic_cluster1, 'Neartic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster1.jpg'), width = get_image_width(nearactic_cluster1), height = get_image_height(nearactic_cluster1), units = "mm")
nearactic_cluster2 = nearctic_cities %>% filter(cluster == 2)
plot_city_cluster(nearactic_cluster2, 'Neartic cluster 2')
Warning: cannot compute exact p-value with tiesWarning: cannot compute exact p-value with tiesWarning: cannot compute exact p-value with tiesWarning: cannot compute exact p-value with ties
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster2.jpg'), width = get_image_width(nearactic_cluster2), height = get_image_height(nearactic_cluster2), units = "mm")
nearactic_cluster3 = nearctic_cities %>% filter(cluster == 3)
plot_city_cluster(nearactic_cluster3, 'Neartic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster3.jpg'), width = get_image_width(nearactic_cluster3), height = get_image_height(nearactic_cluster3), units = "mm")
nearactic_cluster4 = nearctic_cities %>% filter(cluster == 4)
plot_city_cluster(nearactic_cluster4, 'Neartic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neartic_cluster4.jpg'), width = get_image_width(nearactic_cluster4), height = get_image_height(nearactic_cluster4), units = "mm")
neotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Neotropic')
neotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Culiacán" "Guadalajara" "Morelia" "Acapulco" "Querétaro"
[6] "Cuernavaca" "Puebla" "Oaxaca" "Xalapa" "Veracruz"
[11] "Tuxtla Gutiérrez" "Villahermosa" "Guatemala City" "San Salvador" "San Pedro Sula"
[16] "Mérida" "Tegucigalpa" "Managua" "San José" "Cancún"
[21] "Guayaquil" "Chiclayo" "Panama City" "Trujillo" "Quito"
[26] "Havana" "Cali" "Lima" "Pereira" "Miami"
[31] "MedellÃn" "Ibagué" "Cartagena" "Kingston" "Bogota"
[36] "Barranquilla" "Bucaramanga" "Cúcuta" "Maracaibo" "Arequipa"
[41] "Barquisimeto" "Santo Domingo" "Maracay" "El Alto [La Paz]" "Caracas"
[46] "Cochabamba" "Viña del Mar [ValparaÃso]" "RÃo Piedras [San Juan]" "Barcelona" "Concepción"
[51] "Santiago" "Mendoza" "Salta" "Cordoba" "Asuncion"
[56] "Buenos Aires" "La Plata" "Ciudad del Este" "Montevideo" "Mar del Plata"
[61] "Porto Alegre" "São Paulo" "Santos" "Sao Jose dos Campos"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
neotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% neotropic_cities_community_data$city_id))
Run 0 stress 0.134619
Run 1 stress 0.1405831
Run 2 stress 0.1402343
Run 3 stress 0.1346226
... Procrustes: rmse 0.001163646 max resid 0.007319248
... Similar to previous best
Run 4 stress 0.1414222
Run 5 stress 0.134433
... New best solution
... Procrustes: rmse 0.00682457 max resid 0.04571633
Run 6 stress 0.1344509
... Procrustes: rmse 0.001186994 max resid 0.006732139
... Similar to previous best
Run 7 stress 0.134433
... New best solution
... Procrustes: rmse 0.00001887778 max resid 0.00006224176
... Similar to previous best
Run 8 stress 0.1406946
Run 9 stress 0.1346225
... Procrustes: rmse 0.007127633 max resid 0.04569628
Run 10 stress 0.1346406
... Procrustes: rmse 0.007281613 max resid 0.04569614
Run 11 stress 0.134433
... Procrustes: rmse 0.00009355504 max resid 0.0002972054
... Similar to previous best
Run 12 stress 0.1405822
Run 13 stress 0.1405822
Run 14 stress 0.1346226
... Procrustes: rmse 0.007186606 max resid 0.04568917
Run 15 stress 0.1346226
... Procrustes: rmse 0.007196645 max resid 0.04566584
Run 16 stress 0.134619
... Procrustes: rmse 0.0068466 max resid 0.04566077
Run 17 stress 0.134433
... New best solution
... Procrustes: rmse 0.00001770911 max resid 0.00005326707
... Similar to previous best
Run 18 stress 0.1346369
... Procrustes: rmse 0.006925503 max resid 0.04573761
Run 19 stress 0.141222
Run 20 stress 0.1348046
... Procrustes: rmse 0.006551197 max resid 0.04607941
*** Best solution repeated 1 times
neotropic_cities_nmds
scree_plot(neotropic_cities_nmds)
neotropic_cities = cluster_cities(city_nmds = neotropic_cities_nmds, cities_community_data = neotropic_cities_community_data, centers = 5)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(neotropic_cities)
neotropic_biomes = resolve[resolve$REALM == 'Neotropic',c('REALM')]
ggplot() +
geom_sf(data = neotropic_biomes, aes(geometry = geometry)) +
geom_sf(data = neotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_clusters.jpg'))
Saving 7.29 x 4.51 in image
test_value('Neotropic MNTD', neotropic_cities$mntd_normalised)
[1] "Neotropic MNTD mean 0.53 "
test_value('Neotropic beak gape FDiv', neotropic_cities$gape_width_fdiv_normalised)
[1] "Neotropic beak gape FDiv mean 0.44 "
test_value('Neotropic HWI FDiv', neotropic_cities$handwing_index_fdiv_normalised)
[1] "Neotropic HWI FDiv mean 0.53 "
nrow(neotropic_cities)
[1] 64
neotropic_cluster1 = neotropic_cities %>% filter(cluster == 1)
plot_city_cluster(neotropic_cluster1, 'Neotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster1.jpg'), width = get_image_width(neotropic_cluster1), height = get_image_height(neotropic_cluster1), units = "mm")
neotropic_cluster2 = neotropic_cities %>% filter(cluster == 2)
plot_city_cluster(neotropic_cluster2, 'Neotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster2.jpg'), width = get_image_width(neotropic_cluster2), height = get_image_height(neotropic_cluster2), units = "mm")
neotropic_cluster3 = neotropic_cities %>% filter(cluster == 3)
plot_city_cluster(neotropic_cluster3, 'Neotropic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster3.jpg'), width = get_image_width(neotropic_cluster3), height = get_image_height(neotropic_cluster3), units = "mm")
neotropic_cluster4 = neotropic_cities %>% filter(cluster == 4)
plot_city_cluster(neotropic_cluster4, 'Neotropic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster4.jpg'), width = get_image_width(neotropic_cluster4), height = get_image_height(neotropic_cluster4), units = "mm")
neotropic_cluster5 = neotropic_cities %>% filter(cluster == 5)
plot_city_cluster(neotropic_cluster5, 'Neotropic cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'neotropic_cluster5.jpg'), width = get_image_width(neotropic_cluster5), height = get_image_height(neotropic_cluster5), units = "mm")
palearctic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Palearctic')
palearctic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Lisbon" "Porto" "Marrakesh" "Seville" "Dublin" "Málaga"
[7] "Madrid" "Glasgow" "Bilbao" "Liverpool" "Bristol" "Manchester"
[13] "Birmingham" "Leeds" "Newcastle upon Tyne" "Sheffield" "Nottingham" "Valencia"
[19] "London" "Toulouse" "Paris" "Barcelona" "Rotterdam [The Hague]" "Brussels"
[25] "Amsterdam" "Lyon" "Marseille" "Dusseldorf" "Nice" "Frankfurt am Main"
[31] "Zurich" "Oslo" "Stuttgart" "Hamburg" "Genoa" "Nuremberg"
[37] "Copenhagen" "Munich" "Berlin" "Dresden" "Rome" "Prague"
[43] "Stockholm" "Poznan" "Vienna" "Wroclaw" "Zagreb" "Gdansk"
[49] "Budapest" "Krakow" "Warsaw" "Helsinki" "Riga" "Belgrade"
[55] "Lviv" "Sofia" "Thessaloniki" "Saint Petersburg" "Minsk" "Athens"
[61] "Kyiv" "Istanbul" "Odesa" "Samsun" "Luxor" "Tel Aviv"
[67] "Jerusalem" "Tbilisi" "Yerevan" "Kuwait City" "Doha" "Abu Dhabi"
[73] "Dubai" "Bishkek"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
palearctic_cities_nmds = community_nmds(communities %>% filter(city_id %in% palearctic_cities_community_data$city_id))
Run 0 stress 0.04961857
Run 1 stress 0.07367425
Run 2 stress 0.05453968
Run 3 stress 0.0712512
Run 4 stress 0.05826357
Run 5 stress 0.06338732
Run 6 stress 0.06304725
Run 7 stress 0.08618201
Run 8 stress 0.07925636
Run 9 stress 0.05647129
Run 10 stress 0.09129792
Run 11 stress 0.0497659
... Procrustes: rmse 0.005703935 max resid 0.01420466
Run 12 stress 0.05000939
... Procrustes: rmse 0.06743721 max resid 0.2226243
Run 13 stress 0.07603764
Run 14 stress 0.0500095
... Procrustes: rmse 0.06741103 max resid 0.2225553
Run 15 stress 0.096573
Run 16 stress 0.0556821
Run 17 stress 0.05115583
Run 18 stress 0.08625966
Run 19 stress 0.06766821
Run 20 stress 0.06829368
Run 21 stress 0.07723847
Run 22 stress 0.05236354
Run 23 stress 0.05642252
Run 24 stress 0.05235961
Run 25 stress 0.06829346
Run 26 stress 0.05353311
Run 27 stress 0.06421851
Run 28 stress 0.05706493
Run 29 stress 0.0512458
Run 30 stress 0.05479168
*** Best solution was not repeated -- monoMDS stopping criteria:
14: no. of iterations >= maxit
14: stress ratio > sratmax
2: scale factor of the gradient < sfgrmin
palearctic_cities_nmds
scree_plot(palearctic_cities_nmds)
Warning: Quick-TRANSfer stage steps exceeded maximum (= 3700)
palearctic_cities = cluster_cities(city_nmds = palearctic_cities_nmds, cities_community_data = palearctic_cities_community_data, centers = 7)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(palearctic_cities)
palearctic_biomes = st_crop(resolve[resolve$REALM == 'Palearctic',c('REALM')], xmin = -30, ymin = 20, xmax = 80, ymax = 65)
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries
ggplot() +
geom_sf(data = palearctic_biomes, aes(geometry = geometry)) +
geom_sf(data = palearctic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_clusters.jpg'))
Saving 7.29 x 4.51 in image
test_value('Palearctic MNTD', palearctic_cities$mntd_normalised)
[1] "Palearctic MNTD mean 0.58 *"
test_value('Palearctic beak gape FDiv', palearctic_cities$gape_width_fdiv_normalised)
[1] "Palearctic beak gape FDiv mean 0.86 ***"
test_value('Palearctic HWI FDiv', palearctic_cities$handwing_index_fdiv_normalised)
[1] "Palearctic HWI FDiv mean 0.51 "
nrow(palearctic_cities)
[1] 74
palearctic_cluster1 = palearctic_cities %>% filter(cluster == 1)
plot_city_cluster(palearctic_cluster1, 'Palearctic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster1.jpg'), width = get_image_width(palearctic_cluster1), height = get_image_height(palearctic_cluster1), units = "mm")
palearctic_cluster2 = palearctic_cities %>% filter(cluster == 2)
plot_city_cluster(palearctic_cluster2, 'Palearctic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster2.jpg'), width = get_image_width(palearctic_cluster2), height = get_image_height(palearctic_cluster2), units = "mm")
palearctic_cluster3 = palearctic_cities %>% filter(cluster == 3)
plot_city_cluster(palearctic_cluster3, 'Palearctic cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster3.jpg'), width = get_image_width(palearctic_cluster3), height = get_image_height(palearctic_cluster3), units = "mm")
palearctic_cluster4 = palearctic_cities %>% filter(cluster == 4)
plot_city_cluster(palearctic_cluster4, 'Palearctic cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster4.jpg'), width = get_image_width(palearctic_cluster4), height = get_image_height(palearctic_cluster4), units = "mm")
palearctic_cluster5 = palearctic_cities %>% filter(cluster == 5)
nrow(palearctic_cluster5)
[1] 50
plot_city_cluster(palearctic_cluster5[1:25,], 'Palearctic cluster 5a')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster5.jpg'), width = get_image_width(palearctic_cluster5), height = get_image_height(palearctic_cluster5), units = "mm")
plot_city_cluster(palearctic_cluster5[25:50,], 'Palearctic cluster 5b')
Warning: cannot compute exact p-value with ties
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster5b.jpg'), width = get_image_width(palearctic_cluster5), height = get_image_height(palearctic_cluster5), units = "mm")
palearctic_cluster6 = palearctic_cities %>% filter(cluster == 6)
plot_city_cluster(palearctic_cluster6, 'Palearctic cluster 6')
Warning: cannot compute exact p-value with ties
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster6.jpg'), width = get_image_width(palearctic_cluster6), height = get_image_height(palearctic_cluster6), units = "mm")
palearctic_cluster7 = palearctic_cities %>% filter(cluster == 7)
plot_city_cluster(palearctic_cluster7, 'Palearctic cluster 7')
Warning: cannot compute exact p-value with ties
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'palearctic_cluster7.jpg'), width = get_image_width(palearctic_cluster7), height = get_image_height(palearctic_cluster7), units = "mm")
afrotropic_cities_community_data = community_data_metrics %>% filter(core_realm == 'Afrotropic')
afrotropic_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Cape Town" "Johannesburg" "Pretoria" "Kigali" "Kampala" "Arusha" "Nairobi" "Addis Ababa" "Antananarivo"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
afrotropic_cities_nmds = community_nmds(communities %>% filter(city_id %in% afrotropic_cities_community_data$city_id))
Run 0 stress 0.00009014786
Run 1 stress 0.0007443902
Run 2 stress 0.000377487
... Procrustes: rmse 0.002074146 max resid 0.003340817
... Similar to previous best
Run 3 stress 0.00008179516
... New best solution
... Procrustes: rmse 0.0001660244 max resid 0.000361027
... Similar to previous best
Run 4 stress 0.000656661
Run 5 stress 0.3209238
Run 6 stress 0.00009538562
... Procrustes: rmse 0.00003771484 max resid 0.00007382498
... Similar to previous best
Run 7 stress 0.00009646529
... Procrustes: rmse 0.0002653617 max resid 0.0004782319
... Similar to previous best
Run 8 stress 0.00009486991
... Procrustes: rmse 0.0002901629 max resid 0.0004478019
... Similar to previous best
Run 9 stress 0.00009857688
... Procrustes: rmse 0.0004629411 max resid 0.0009418255
... Similar to previous best
Run 10 stress 0.00008547616
... Procrustes: rmse 0.00001297639 max resid 0.00002426197
... Similar to previous best
Run 11 stress 0.00009878898
... Procrustes: rmse 0.0004635066 max resid 0.0009419129
... Similar to previous best
Run 12 stress 0.0005427305
... Procrustes: rmse 0.002940686 max resid 0.004342951
... Similar to previous best
Run 13 stress 0.001227374
Run 14 stress 0.00009793901
... Procrustes: rmse 0.0001331878 max resid 0.0002796824
... Similar to previous best
Run 15 stress 0.00009902405
... Procrustes: rmse 0.0001283685 max resid 0.0002454039
... Similar to previous best
Run 16 stress 0.00009120045
... Procrustes: rmse 0.00002743712 max resid 0.00005350613
... Similar to previous best
Run 17 stress 0.00009094513
... Procrustes: rmse 0.0002505202 max resid 0.00046142
... Similar to previous best
Run 18 stress 0.00009137309
... Procrustes: rmse 0.000123082 max resid 0.0002401071
... Similar to previous best
Run 19 stress 0.00009536048
... Procrustes: rmse 0.00003982921 max resid 0.00006129269
... Similar to previous best
Run 20 stress 0.00006560657
... New best solution
... Procrustes: rmse 0.0001813766 max resid 0.0003094943
... Similar to previous best
*** Best solution repeated 1 times
Warning: stress is (nearly) zero: you may have insufficient data
afrotropic_cities_nmds
scree_plot(afrotropic_cities_nmds)
afrotropic_cities = cluster_cities(city_nmds = afrotropic_cities_nmds, cities_community_data = afrotropic_cities_community_data, centers = 2)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(afrotropic_cities)
afrotropic_biomes = resolve[resolve$REALM == 'Afrotropic',c('REALM')]
ggplot() +
geom_sf(data = afrotropic_biomes, aes(geometry = geometry)) +
geom_sf(data = afrotropic_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_clusters.jpg'))
Saving 7.29 x 4.51 in image
test_value('Afrotropic MNTD', afrotropic_cities$mntd_normalised)
[1] "Afrotropic MNTD mean 0.14 *"
test_value('Afrotropic beak gape FDiv', afrotropic_cities$gape_width_fdiv_normalised)
[1] "Afrotropic beak gape FDiv mean 0.41 "
test_value('Afrotropic HWI FDiv', afrotropic_cities$handwing_index_fdiv_normalised)
[1] "Afrotropic HWI FDiv mean 0.61 "
nrow(afrotropic_cities)
[1] 9
afrotropic_cluster1 = afrotropic_cities %>% filter(cluster == 1)
plot_city_cluster(afrotropic_cluster1, 'Afrotropic cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster1.jpg'), width = get_image_width(afrotropic_cluster1), height = get_image_height(afrotropic_cluster1), units = "mm")
afrotropic_cluster2 = afrotropic_cities %>% filter(cluster == 2)
plot_city_cluster(afrotropic_cluster2, 'Afrotropic cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'afrotropic_cluster2.jpg'), width = get_image_width(afrotropic_cluster2), height = get_image_height(afrotropic_cluster2), units = "mm")
indomalayan_cities_community_data = community_data_metrics %>% filter(core_realm == 'Indomalayan')
indomalayan_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Srinagar" "Jamnagar" "Jammu" "Rajkot" "Bikaner" "Jodhpur" "Jalandhar"
[8] "Ahmedabad" "Bhavnagar" "Ludhiana" "Anand" "Udaipur" "Surat" "Vadodara"
[15] "Ajmer" "Chandigarh" "Vasai-Virar" "Mumbai" "Jaipur" "Delhi [New Delhi]" "Nashik"
[22] "Dehradun" "Kota" "Pune" "Haridwar" "Dhule" "Ujjain" "Indore"
[29] "Ahmadnagar" "Kolhapur" "Jalgaon" "Agra" "Aurangabad" "Sangli" "Belagavi"
[36] "Gwalior" "Budaun" "Bareilly" "Dharwad" "Bhopal" "Bhind" "Mangaluru"
[43] "Solapur" "Vijayapura" "Akola" "Latur" "Kannur" "Davanagere" "Thalassery"
[50] "Amravati" "Kalaburagi" "Kozhikode" "Guruvayur" "Malappuram" "Lucknow" "Thrissur"
[57] "Mysuru" "Kochi" "Alappuzha" "Nagpur" "Kollam" "Jabalpur" "Ettumanoor"
[64] "Hyderabad" "Coimbatore" "Bengaluru" "Thiruvananthapuram" "Tiruppur" "Faizabad" "Erode"
[71] "Prayagraj" "Pratapgarh" "Salem" "Dindigul" "Madurai" "Tiruchirappalli" "Durg"
[78] "Vellore" "Tirupati" "Raipur" "Bilaspur" "Vijayawada" "Puducherry" "Chennai"
[85] "Kathmandu" "Colombo" "Rajamahendravaram" "Patna" "Kandy" "Bihar Sharif" "Visakhapatnam"
[92] "Ranchi" "Brahmapur" "Jamshedpur" "Darjeeling" "Siliguri" "Cuttack" "Bhubaneshwar"
[99] "Jalpaiguri" "Berhampore" "Kolkata" "Krishnanagar" "Guwahati [Dispur]" "Agartala" "Silchar"
[106] "Dimapur" "Bangkok" "George Town" "Kuala Lumpur" "Phnom Penh" "Singapore" "Hong Kong"
[113] "Sha Tin" "Hsinchu" "Taichung" "New Taipei [Taipei]" "Tainan" "Denpasar" "Kaohsiung"
[120] "Kota Kinabalu"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
indomalayan_cities_nmds = community_nmds(communities %>% filter(city_id %in% indomalayan_cities_community_data$city_id))
Run 0 stress 0.1190668
Run 1 stress 0.1224401
Run 2 stress 0.1394961
Run 3 stress 0.1161598
... New best solution
... Procrustes: rmse 0.02798297 max resid 0.2346031
Run 4 stress 0.1540252
Run 5 stress 0.1175538
Run 6 stress 0.1167657
Run 7 stress 0.1384412
Run 8 stress 0.1241275
Run 9 stress 0.1163498
... Procrustes: rmse 0.025336 max resid 0.2481842
Run 10 stress 0.1153501
... New best solution
... Procrustes: rmse 0.008287738 max resid 0.08217213
Run 11 stress 0.1199235
Run 12 stress 0.1172201
Run 13 stress 0.151756
Run 14 stress 0.1246307
Run 15 stress 0.134023
Run 16 stress 0.1580175
Run 17 stress 0.1229198
Run 18 stress 0.1502633
Run 19 stress 0.1189366
Run 20 stress 0.117072
Run 21 stress 0.1458634
Run 22 stress 0.137476
Run 23 stress 0.1186188
Run 24 stress 0.1637221
Run 25 stress 0.1191928
Run 26 stress 0.1164559
Run 27 stress 0.142366
Run 28 stress 0.1379337
Run 29 stress 0.1191853
Run 30 stress 0.1266232
*** Best solution was not repeated -- monoMDS stopping criteria:
29: stress ratio > sratmax
1: scale factor of the gradient < sfgrmin
indomalayan_cities_nmds
scree_plot(indomalayan_cities_nmds)
indomalayan_cities = cluster_cities(city_nmds = indomalayan_cities_nmds, cities_community_data = indomalayan_cities_community_data, centers = 5)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(indomalayan_cities)
indomalayan_biomes = resolve[resolve$REALM == 'Indomalayan',c('REALM')]
ggplot() +
geom_sf(data = indomalayan_biomes, aes(geometry = geometry)) +
geom_sf(data = indomalayan_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_clusters.jpg'))
Saving 7.29 x 4.51 in image
test_value('Indomalayan MNTD', indomalayan_cities$mntd_normalised)
[1] "Indomalayan MNTD mean 0.44 ***"
test_value('Indomalayan beak gape FDiv', indomalayan_cities$gape_width_fdiv_normalised)
[1] "Indomalayan beak gape FDiv mean 0.52 "
test_value('Indomalayan HWI FDiv', indomalayan_cities$handwing_index_fdiv_normalised)
[1] "Indomalayan HWI FDiv mean 0.82 ***"
nrow(indomalayan_cities)
[1] 120
indomalayan_cluster1 = indomalayan_cities %>% filter(cluster == 1)
plot_city_cluster(indomalayan_cluster1, 'Indomalayan cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster1.jpg'), width = get_image_width(indomalayan_cluster1), height = get_image_height(indomalayan_cluster1), units = "mm")
indomalayan_cluster2 = indomalayan_cities %>% filter(cluster == 2)
plot_city_cluster(indomalayan_cluster2, 'Indomalayan cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster2.jpg'), width = get_image_width(indomalayan_cluster2), height = get_image_height(indomalayan_cluster2), units = "mm")
indomalayan_cluster3 = indomalayan_cities %>% filter(cluster == 3)
plot_city_cluster(indomalayan_cluster3, 'Indomalayan cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster3.jpg'), width = get_image_width(indomalayan_cluster3), height = get_image_height(indomalayan_cluster3), units = "mm")
indomalayan_cluster4 = indomalayan_cities %>% filter(cluster == 4)
plot_city_cluster(indomalayan_cluster4, 'Indomalayan cluster 4')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster4.jpg'), width = get_image_width(indomalayan_cluster4), height = get_image_height(indomalayan_cluster4), units = "mm")
indomalayan_cluster5 = indomalayan_cities %>% filter(cluster == 5)
plot_city_cluster(indomalayan_cluster5, 'Indomalayan cluster 5')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'indomalayan_cluster5.jpg'), width = get_image_width(indomalayan_cluster5), height = get_image_height(indomalayan_cluster5), units = "mm")
australasia_cities_community_data = community_data_metrics %>% filter(core_realm == 'Australasia')
australasia_cities_community_data %>% dplyr::select(city_name) %>% distinct() %>% as.list()
$city_name
[1] "Perth" "Adelaide" "Melbourne" "Sydney" "Brisbane" "Auckland"
attr(,"na.action")
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
1 56 81 87 90 92 94 96 98 100 102 103 104 106 112 115 116 117 118 119 121 215
attr(,"class")
[1] "omit"
australasia_cities_nmds = community_nmds(communities %>% filter(city_id %in% australasia_cities_community_data$city_id))
Run 0 stress 0
Run 1 stress 0
... Procrustes: rmse 0.1120339 max resid 0.1622729
Run 2 stress 0.1572548
Run 3 stress 0
... Procrustes: rmse 0.07245479 max resid 0.1129183
Run 4 stress 0
... Procrustes: rmse 0.04603582 max resid 0.07090683
Run 5 stress 0
... Procrustes: rmse 0.2182675 max resid 0.295989
Run 6 stress 0
... Procrustes: rmse 0.133075 max resid 0.2072423
Run 7 stress 0.0000995744
... Procrustes: rmse 0.1149167 max resid 0.1742941
Run 8 stress 0
... Procrustes: rmse 0.1694089 max resid 0.2635754
Run 9 stress 0
... Procrustes: rmse 0.1063225 max resid 0.1577545
Run 10 stress 0
... Procrustes: rmse 0.1313955 max resid 0.2135724
Run 11 stress 0
... Procrustes: rmse 0.1251736 max resid 0.1780552
Run 12 stress 0
... Procrustes: rmse 0.05487419 max resid 0.08977988
Run 13 stress 0.1572548
Run 14 stress 0.00006459844
... Procrustes: rmse 0.1484605 max resid 0.2073383
Run 15 stress 0.00008456015
... Procrustes: rmse 0.1937182 max resid 0.3430883
Run 16 stress 0
... Procrustes: rmse 0.09381547 max resid 0.1446598
Run 17 stress 0.1572548
Run 18 stress 0.00001782659
... Procrustes: rmse 0.132877 max resid 0.2076538
Run 19 stress 0
... Procrustes: rmse 0.115211 max resid 0.1946323
Run 20 stress 0
... Procrustes: rmse 0.1542908 max resid 0.2134235
Run 21 stress 0.1572548
Run 22 stress 0
... Procrustes: rmse 0.1239607 max resid 0.190374
Run 23 stress 0.00007261303
... Procrustes: rmse 0.123741 max resid 0.1717834
Run 24 stress 0.00006062953
... Procrustes: rmse 0.1354649 max resid 0.2269844
Run 25 stress 0.00008963805
... Procrustes: rmse 0.1180163 max resid 0.1705588
Run 26 stress 0.000009458018
... Procrustes: rmse 0.1301844 max resid 0.1953992
Run 27 stress 0.2349416
Run 28 stress 0.0000881637
... Procrustes: rmse 0.07152338 max resid 0.1085942
Run 29 stress 0
... Procrustes: rmse 0.09917667 max resid 0.1266893
Run 30 stress 0
... Procrustes: rmse 0.08881426 max resid 0.1290499
*** Best solution was not repeated -- monoMDS stopping criteria:
25: stress < smin
1: stress ratio > sratmax
4: scale factor of the gradient < sfgrmin
Warning: stress is (nearly) zero: you may have insufficient data
australasia_cities_nmds
scree_plot(australasia_cities_nmds)
australasia_cities = cluster_cities(city_nmds = australasia_cities_nmds, cities_community_data = australasia_cities_community_data, centers = 3)
Joining with `by = join_by(city_id)`
plot_nmds_clusters(australasia_cities)
australasia_biomes = st_crop(resolve[resolve$REALM == 'Australasia',c('REALM')], xmin = 80, ymin = -70, xmax = 180, ymax = 20)
although coordinates are longitude/latitude, st_intersection assumes that they are planar
Warning: attribute variables are assumed to be spatially constant throughout all geometries
ggplot() +
geom_sf(data = australasia_biomes, aes(geometry = geometry)) +
geom_sf(data = australasia_cities, aes(geometry = geometry, color = cluster))
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'australasia_clusters.jpg'))
Saving 7.29 x 4.51 in image
test_value('Australasia MNTD', australasia_cities$mntd_normalised)
[1] "Australasia MNTD mean 0.48 "
test_value('Australasia beak gape FDiv', australasia_cities$gape_width_fdiv_normalised)
[1] "Australasia beak gape FDiv mean 0.5 "
test_value('Australasia HWI FDiv', australasia_cities$handwing_index_fdiv_normalised)
[1] "Australasia HWI FDiv mean 0.82 "
nrow(australasia_cities)
[1] 6
australasia_cluster1 = australasia_cities %>% filter(cluster == 1)
plot_city_cluster(australasia_cluster1, 'Australasia cluster 1')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'australasia_cluster1.jpg'), width = get_image_width(australasia_cluster1), height = get_image_height(australasia_cluster1), units = "mm")
australasia_cluster2 = australasia_cities %>% filter(cluster == 2)
plot_city_cluster(australasia_cluster2, 'Australasia cluster 2')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'australasia_cluster2.jpg'), width = get_image_width(australasia_cluster2), height = get_image_height(australasia_cluster2), units = "mm")
australasia_cluster3 = australasia_cities %>% filter(cluster == 3)
plot_city_cluster(australasia_cluster3, 'Australasia cluster 3')
ggsave(filename(REGION_DEEP_DIVE_FIGURES_OUTPUT, 'australasia_cluster3.jpg'), width = get_image_width(australasia_cluster3), height = get_image_height(australasia_cluster3), units = "mm")